Cluster approximations for infection dynamics on random networks 
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In this paper, we consider a simple stochastic epidemic model on large regular random graphs and 
the stochastic process that corresponds to this dynamics in the standard pair approximation. Using 
the fact that the nodes of a pair are unlikely to share neighbors, we derive the master equation for 
this process and obtain from the system size expansion the power spectrum of the fluctuations in 
the quasi-stationary state. We show that whenever the pair approximation deterministic equations 
give an accurate description of the behavior of the system in the thermodynamic limit, the power 
spectrum of the fluctuations measured in long simulations is well approximated by the analytical 
power spectrum. If this assumption breaks down, then the cluster approximation must be carried 
out beyond the level of pairs. We construct an uncorrelated triplet approximation that captures the 
behavior of the system in a region of parameter space where the pair approximation fails to give a 
good quantitative or even qualitative agreement. For these parameter values, the power spectrum 
of the fluctuations in finite systems can be computed analytically from the master equation of the 
corresponding stochastic process. 

PACS numbers: 87.10.Mn, 05.10.Gg, 87.10.Ca, 87.10.Rt 



I. INTRODUCTION 

Stochastic models on lattices are an old subject in 
statistical physics, and cluster mean-field theories have 
been developed and used both in the context of equi- 
librium and non-equilibrium problems [l|, Q ■ These are 
mean-field approximations that are formulated in terms 
of the n-site joint probabilities. Since for any fixed n 
the evolution equations for these probabilities do not in 
general form a closed set, the (n + l)-site probabilities 
are expressed in terms of the lower order probabilities 
according to some closure assumption. At the lowest 
nontrivial order of truncation, this method corresponds 
to the well known pair approximation (PA) according 
to which the triplet probabilities p(a, b, c) are factorizcd 
as p(a,b,c) = p{a,b)p(c\b), where a, 6, c denote lattice 
nodes' states and p(c\b) is the conditional probability of 
having state c in neighbor of a node in state b. The PA is 
exact on the Bethe lattice or Cayley tree, a mathematical 
model which cannot be realized as a physical system or 
simulated computationally. In finite lattices, in particu- 
lar in the d-dimensional regular lattice with first neighbor 
interactions often considered in these models, the PA is 
used as the simplest analytical description that includes 
an explicit representation of spatial correlations. Despite 
being quantitatively inaccurate, it provides insight and 
qualitative information about the system. 

More recently, the interest in this approximation pro- 
cedure shifted from the exploration of conceptual models 
to its applications in several problems of population dy- 
namics. The spread of a virus is an example of a dynamic 
process occurring on a discrete spatial arrangement that 
was modeled in the traditional literature with the 1-site 
or mean- field approximation (MFA) 0, 0|. While the 
MFA reasonably reproduces the spreading behavior for 



topologies where the number of connections per node 
is either high or strongly fluctuating and for those that 
show small- world features, it is highly inaccurate for lat- 
tice and in general network structured populations. The 
PA has become very popular in lattice-based stochastic 
models of ecological, epidemic and evolutionary game dy- 
namics 0, 0, H, H, 2, 11, H, HI 14, 15 1, and several 
modifications of the PA have been proposed for particular 
graphs and dynamic rules that lead to better quantitative 
agreement [iB, E3, H E 0, ED, 0] • Some of these 
modifications are based on different closure assumptions 
at the level of pairs, while others depend on including 
higher order clusters. 

By contrast, the PA is expected to perform well, even 
quantitatively, for stochastic models on large regular ran- 
dom graphs (RRGs) which are random networks of fixed 
connectivity per node, or degree [24[ . This is because of 
the RRG's statistical properties, namely the short loop 
density of the graph tending to zero in the limit of large 
graph size [25(, so that a RRG may be seen locally as a 
Bethe lattice of the same degree. 

In this paper, we consider the susceptible-infcctive- 
recovered-susceptible (SIRS) stochastic epidemic model 
on simple RRGs, that is RRGs with no loops formed 
by one or two edges. We construct a detailed stochastic 
model based on the PA that captures the behavior of this 
dynamics in finite systems, in the sense that the power 
spectrum of the fluctuations computed analytically from 
the model matches the numerical power spectrum mea- 
sured along the simulation runs whenever the averaged 
dynamics of the densities is well approximated by the so- 
lutions of the PA equations in the thermodynamic limit. 
This happens in a large region of parameter space be- 
cause the quality of the approximation becomes poor only 
when the recovery rate is much larger than the rate of 
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loss of immunity. The coupling between dynamics and 
graph structure intervenes in the microscopic description 
of the system through the different transitions that it is 
necessary to consider and the corresponding transition 
rates. 

When recovery is much faster than loss of immunity 
and the PA fails in the thermodynamics limit, the clus- 
ter approximation must be carried out beyond the level of 
pairs. We show that a triplet approximation (TA) with a 
standard closure assumption is adequate to describe the 
behavior of the system up to recovery rates two orders of 
magnitude larger than the immunity waning rates. This 
parameter range is of interest in applications to childhood 
infectious diseases [3, 0] • In order to describe the behav- 
ior of finite systems in this region, a detailed stochastic 
model suitable for RRGs can be built on the basis of the 
TA following the procedure described for the PA. 



II. DETERMINISTIC AND STOCHASTIC 
FRAMEWORKS IN THE PAIR 
APPROXIMATION 

In this section we propose a stochastic model of the 
SIRS epidemic process in the PA. The effects of the inclu- 
sion of an implicit representation of spatial dependence 
in a stochastic model have been recently studied [27( . An 
extension of the analysis of stochastic fluctuations from 
non-spatial models to the case of models on regular struc- 
tures such as d-dimensional hypercubic lattices has also 
been performed 28] . The present study elaborates on the 
stochastic models with implicit spatial dependence by in- 
cluding a detailed microscopic description of the transi- 
tions between the states of the nodes and the states of 
the pairs of the nearest neighbors, that can be treated 
analytically. The results obtained both for the averaged 
dynamics that describes the behavior of the system in 
the thermodynamic limit and for the spectrum of the 
fluctuations in the quasi-stationary state of finite systems 
agree well with the results of Monte Carlo simulations on 
RRGs in a large parameter range. Technically, the PA 
approach is quite similar to that previously developed for 
well-mixed systems also known as the MFA approach [2§| 
while having substantial differences which can be easily 
understood as soon as we review a few facts. 

The full set of deterministic equations describing the 
SIRS process in the MFA is written in terms of the den- 
sities of susceptible and infected populations (we will 
equivalently use the term "probability" in what follows). 
These differential equations are deduced on the assump- 
tion of uncorrelated nodes in the limit of infinite sys- 
tem size and constitute an approximation for the tran- 
sient and quasi-stationary behaviors of large spatially ex- 
tended systems. Accordingly, the associated stochastic 
model has two classes of individuals, infected and sus- 
ceptible, as its independent dynamical variables, and a 



particular realization of the model at a given time t con- 
sists of mi susceptible individuals, mi infected individu- 
als and (N — m± — mi) recovered individuals, where N 
is the total population size or the number of nodes in a 
graph. 

The deterministic formulation of the SIRS model in the 
PA couples the dynamics of the node and pair densities in 
the thermodynamic limit. The purpose of the inclusion 
of pair densities is to improve the level of approximation 
in the description of spatially explicit simulations of the 
stochastic model. Consider the stochastic SIRS process 
on a RRG with fixed degree per node k (RRG-fc) and N 
nodes. After the initial distribution of the nodes among 
the classes of susceptible (£>), infected (I) and recovered 
(R) individuals, let the state of the system evolve in time 
according to asynchronous update of the events of infec- 
tion, recovery and immunity waning. Namely, infected 
individuals recover at rate S [I — ► R] , immunity of recov- 
ered individuals ceases at rate 7 [R — > S], and infection in 
the susceptible-infected pairs of the nearest neighbors oc- 
curs with rate A [SI — > II]. Each executed event taken 
separately changes the state of only one node, has its 
side-effects on the state of the pairs the node forms with 
its neighbors, propagates to larger configurations and fi- 
nally to the whole graph. Consequently, frequency distri- 
butions or numbers of various correlations present in the 
graph, among which pairs are the simplest ones, change. 
To construct a stochastic model in the PA the transition 
probabilities of different events have to be approximated 
purely in terms of pairs. 

More specifically, to represent independent classes of 
the constituents of the stochastic system we can choose 
five independent numbers subject to the constraints on 
the number of nodes and pairs in the graph. These can 
be five pair variables alone, one node and four pair vari- 
ables or two node and three pair variables. We shall 
adopt the last variant so that taking the limit N — ► 00 
in the stochastic PA model we recover the PA-SIRS de- 
terministic equations in the form given in Ref. [27| . We 
introduce the following notation for the independent vari- 
ables: mi and m^ arc the numbers of susceptible and 
of infected nodes and 7713, 7774 and 777,5 ar e the num- 
bers of susceptible-infected, susceptible-recovered and 
recovered-infected pairs of the nearest neighbor nodes, 
respectively, at time t. The full set of the five variables 
is denoted shortly as m = (mi, TO2, 7713, 777,4, 7775). The 
remaining variables are not independent. They can be 
computed from the constraints as follows: 



N = y^m a , 

a 

km a = map + 2m aa , 



(1) 

(2) 



where m a is the number of nodes in state a, m a p = mp a 
is the number of pairs of the nearest neighbor nodes in 
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FIG. 1: Schematic representation of a node in state a and its 
nearest neighbors in states pi, where i = 1, . . . , 4 and a, Pi £ 
{S,I,R}, in a RRG-4. 
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FIG. 2: Example of a configuration before and after recovery 
of the central node. The total number of recovery configura- 
tions is 15. 



states a, (3 € {S, I, R}. In this notation, mi and m 2 
equal ms and mj, and 777,3, 7714 and 7715 equal mg/, itisr 
and mjj;, respectively. Note that in Eq. the factor 
of 2 comes from the fact that aa pairs must be counted 
twice. Further, the number of constituents in the classes, 
which are now interpreted as classes of individuals in 
specific states and classes of fixed contacts established 
amongst individuals in specific states, evolves according 
to the rule: whichever event is executed, a transition of 
one individual between the classes of individuals occurs 
simultaneously with transitions of the individual's k con- 
tacts between the classes of contacts. This rule is nothing 
else than the statement of the local modifications on a 
RRG-fc caused by a single event: a change in state of one 
node induces the simultaneous change in states of the k 
pairs whose common member is the node undergoing the 
transition. The coarse grained description where the ef- 
fect of the change in state of a given node on the k pairs 
that it forms is averaged over each pair type has been 
given in Ref. [27|. In this study, we consider a detailed 
stochastic model in which the k pairs are, in general, 
different. 

Next, we calculate the transition rates for the PA-SIRS 
stochastic process corresponding to the full microscopic 
description and yielding the PA deterministic equations 
as the equations of motion for the infinite system size. 
Consider a recovery event [/ — > R] on a RRG-4 as an 
example. During a simulation of the stochastic process 
the rate of this event is equal to the constant rate of re- 
covery, 6, multiplied by the number of infected nodes at 
given time, 7712. As we keep track of the changes in states 
of the pairs formed by the central node that switches 
from infected to recovered, we have to calculate how this 
net rate is distributed among all possible configurations 
the infected node and its first neighbors might form. The 
configurations with a central node in state / are accepted 
as different if they have different number of pairs mj a , 
where a € {S, /, R}, irrespective of their spatial arrange- 
ment. The total number of such configurations for recov- 
ery equals ( n+ ^.~ 1 ), where n — 3 is the number of states 
and k = 4 is coordination number of the graph. 

As it is well known in large sparse RRGs the members 
of a pair are unlikely to share neighbors [25|, [3(| • Having 
no nodes in common suggests that in the first approxi- 



mation the pairs can be considered as independently dis- 
tributed. In mathematical formulation it means that the 
probabilities of having particular configurations follow a 
multinomial law [TtJ ■ To be as general as possible we give 
the conditional probability of a configuration, shown in 
Fig. [TJ by the formula: 

( ^ 

V Ps 

= M f [ Q(m m ;\ (3) 



m a (3 



Here a, fa £ {S, I, R} denote states of the five nodes and 
i = 1, . . . , k. The number of neighbors to be indepen- 
dently distributed is k = 4. Q(fa\a) is the conditional 
probability that given a node in state a its nearest neigh- 
bor is in state fa: 
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Q(fa\a) 
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if a = fa 



(4) 



Summation over the states of all pairs present in a given 
configuration equals the number of the nearest neighbors: 

E m aPi = k - 
aPi 

If there are chains of connections between outer nodes 
of a configuration, especially short ones, the described 
method introduces an error because the pairs can no 
longer be considered independent. This is a crucial point 
that explains the failure of the standard PA for regular 
structures characterized by a huge number of loops of all 
lengths. For sparse RRGs with locally tree-like structure 
we assume the pairs to be uncorrelated, so that the prob- 
ability distribution of a central node's neighbors can still 
be estimated by the multinomial Eq. ([3]). 

To demonstrate an application of Eq. ([3]) we choose a 
configuration before recovery of the central node shown 
in Fig. [2] The probability of this configuration is given 
by: 
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Using Eq. ([2]) and the notation introduced in the begin- 
ning of this section one obtains: 



4(km 2 — ins — ms) ms 




[km 2 f 



According to Fig. [TJ we introduce the total variation 
of the number of nodes in state a and of the number of 
pairs in state afli in a graph as the difference between 
the respective numbers after and before the central node 
switches from state a to another state: 
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In Fig. [21 after recovery of the central node the con- 
figuration changes to that with a recovered node in the 
center inducing transitions in pairs of the nearest neigh- 
bors. The corresponding number variations for the inde- 
pendent variables are Ami — 0, Am 2 = — 1, A7713 = — 1, 
Ami — 1, and Am§ = 3. Finally, we can give the tran- 
sition rate for a recovery event in the center of the con- 
figuration shown in Fig. [5] as follows: 



^-mi,m2-l,m3-l,"i4+l ; ras+3 _ 
mi,m2,m3,m4,Tfi5 '' '-- 



4(fcm2 — TO3 — TO5) 3 77l3 

(fcm 2 ) 4 



where the subscript and the superscript of T denote the 
initial and the final states of the system, respectively. 

The other transition rates for events of recovery and 
immunity waning are straightforward modifications of 

the above. As for infection process [SI — > II] it is consid- 
ered as inherent to a pair of nodes because transmission 
of infection requires a contact between two individuals, 
susceptible and infected. A transition of S node in SI 
pair from susceptible to infected changes both the state 
of this pair and the states of the other three pairs join- 
ing the S node with its three nearest neighbors which 
we assume to be independently distributed in the PA. 
Thus, Eq. ([3]) is still applicable with a = S, $4 — I and 
i = 1 k — 1: 
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where we have used the closure assumption in the stan- 
dard PA, Q(Pi\SI) « Q(fli\S). Now, for instance, the 
approximate transition rate for infection within an SI 
pair of the configuration depicted in Fig. [3] equals: 
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FIG. 3: Example of a configuration before and after infection. 
The total number of infection configurations is 10. 



Having calculated all transition rates as described, we 
are ready to write the master equation for the PA-SIRS 
stochastic process 3l|, |32[ : 



dP(m,t) 
Jt 



E 



\T™,V{m',t)-T™'v{m,t) 



(6) 



where T™, denotes transition rates from other states m! 
to state m and vice versa for T™ . The complete solu- 
tion of this differential-difference equation V(m,t) gives 
the probability of finding the system in state m for all 
allowed sets of integers m„ where i — 1, . . . , 5, at time 
t > subject to the initial, normalization and boundary 
conditions. In general, it is not easy to solve this equa- 
tion analytically but it is quite straightforward to analyze 
it for large but finite N using van Kampen's system size 
expansion [31]. In that spirit, we set: 



mi(t) = NP(S)(t) 
m 2 (t) = NP(I)(t)- 



N Xl (t) , 



(7) 



In both equations, the first macroscopic terms scale 
with the system size N. The functions P(S)(t) = 
lim m\(t)/N and P(I)(t) = lim rri2(t)/N are densi- 

N—>oo N^oc 

ties of susceptible and infected populations which have 
to be adjusted so as to satisfy the deterministic equa- 
tions of motion in the PA. x\{t) and x 2 {i) are the new 
variables which denote stochastic fluctuations around the 
corresponding solutions of the PA deterministic equa- 
tions and replace m\(t) and m 2 (t) 7 respectively. The 
time-dependent transformations [Eq. iff]) ] from mi(t), 
m 2 (t) to Xi(t), x 2 (t) involving functions P(S)(t), P(I)(t) 
come from the fact that one expects, with respect to 
the node variables, the probability distribution function 
V(m, t) to have a sharp peak around the macroscopic val- 
ues mi(t) = NP(S)(t), m 2 {t) = NP{I){t) with a width 
of order of y/N, so that the functions P(S)(t), P(I)(t) 
follow the motion of the peak in time. Wheareas the 
system involves nodes and pairs as constituent elements 
we should carefully define what is meant by stochastic 
fluctuations around the solutions of the PA deterministic 
equations regarding the pair variables. For the fluctua- 
tions of the pair densities we set: 

m 3 (t) = NkP(SI)(t) - 
mA(t) = NkP{SR){t) 
m 5 (t) = NkP(RI)(t) - 



Nkx 3 (t) , 
Nkx 4 (t) , (8) 
Nkx 5 (t) . 
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In the above transformations, the first macroscopic terms 
are of order Nk, and the ansatz for the fluctuations scal- 
ing with the system size is V^/Vfc for fixed parameter fc 
[actually, the scaling of the macroscopic terms in Eqs. 
and © can be found from Eqs. Q and ©]. The 
densities of susceptible-infected, susceptible-recovered 
and recovered-infected pairs are defined as P(SI)(t) = 
lim m 3 (t)/(Nk), P(SR)(t) = lim rndt)/(Nk) and 

N— too N— too 

P(RI)(t) = lim m 5 (t)/(Nk), respectively. 

N— too 

The large- N expansion is carried out using integer step 
operators e.;, where i = 1, . . . , 5, that can be expressed as 
Taylor series involving partial derivatives with respect to 
the node and pair fluctuation variables [3lj [: 
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With the aid of Eqs. ©-©I, the master equation 
be written in the following form: 



(9) 



can 



dH(x,t) 
dt 



E (n^-iW 11 ^)' ( 10 ) 



where the probability distribution function H(x, t) 
V{m, t) and rij are integer powers such that: 



n 



The leading terms in the power series expansion of Eq. 
(fTTJ)) are of order y/N. Collecting and setting them to zero 
yields a set of five first order differential equations. These 
are the PA-SIRS deterministic equations, see the Ap- 
pendix for their explicit form. In next-to-leading-order, 
setting to zero the terms of order N°, one gets a linear 
multivariate Fokker-Planck equation for the probability 
distribution function H(x,t) [3ll. l32j: 



dU 
~dt 



■5> 



d(xjIL) 

dxi 



d 2 n 

13 dxidXn 



(11) 



Here x = [x\, X2, %a, x<i, x§) are stochastic fluctuations of 
the node and pair densities about their endemic steady 
state values in the PA. A is the Jacobian matrix of the 
PA equations linearized about the endemic equilibrium 
solution, see formula (]28[) in the Appendix. B is sym- 
metric internal noise cross correlation matrix derived di- 
rectly from the expansion. For instance, Bn — —SP(I), 
where P(I) stands for the endemic equilibrium value of 
the density of infected individuals in the PA. The differ- 
ence of the detailed stochastic description from the coarse 
grained description considered in Ref. 27] is reflected in 



Since the Fokker-Planck equation is linear, its solution 
II(a;, i) is a multivariate Gaussian distribution completely 
determined by the first and the second moments. How- 
ever, to analyze the fluctuations it is convenient to use 
the equivalent linear multivariate Langevin equation for 
Xi{t) 0,113: 



x i {t)=Y,A ij x j {t)+L i {t) , i,j = 1 5 . (12) 

3 

Li(t) are white random noise terms with the following 
properties: 



(Li®) = , 

(L i (t)L j (t'))=B ij 6{t-t') 



(13) 



The structure of the noise Xi(t) as function of angular 
frequency uj is found from the power spectrum of the 
normalized fluctuations (PSNF) denoted as the averaged 
squared modulus of the Fourier transform of Xi(t): 



where 



Xi(uj) 



/2tt 



-iujt 



(14) 



(15) 



Solving for the Fourier transforms from the linear Eq. 
(|12p and using the correlation function in the frequency 
domain (Li(u>)Lj(u>')) = BijS(uj + u/), the approximate 
analytical expression for the PSNFs about the endemic 
equilibrium solution of the PA equations becomes: 



P i (u;) = J2Mr k 1 (u)B kj Mr i \- 

3,k 



(16) 



the elements of the matrix Bij, where i,j = 3, 4, 5. 



where Mij(oj) — iui 5ij — Aij. The PSNF of the suscepti- 
bles (of the infectives) is then obtained by setting i = 1 
(i = 2, respectively). In this case, the PSNFs are of the 
form p(u>)/q(w), where p(u>) and q(u>) are polynomials in 
u> of order 8 and 10, respectively. Note that the same for- 
mula [Eq. (HU)] is also valid for the PSNFs in the MFA 
taking A as the Jacobian of the MFA-SIRS differential 
equations linearized about the mean-field endemic equi- 
librium solution and the noise cross correlation matrix B 
computed directly from the expansion. 

Figs. S] and [5] compare results of the theory developed 
so far with data of the SIRS stochastic model obtained 
from Monte Carlo simulations on a RRG-4. The RRGs 
are generated using a quick algorithm introduced in Ref. 
[33j ]. The algorithm guarantees that for small degrees 
at least, which is the case here, the RRG-fc on N nodes 
is generated uniformly at random, in the sense that 
all RRG-fc on TV nodes have asymptotically the same 
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probability as N — > oo. In the stochastic simulations, 
the system is set in a random initial condition with 
fixed node and pair densities, after which the states 
of the nodes are updated asynchronously according to 
the events of infection, recovery and loss of immunity. 
Figure H] shows averaged or global behavior of the 
densities as obtained from many runs of numerical 
simulations, namely for each set of given parameter 
values and initial conditions the simulations are averaged 
over 10 3 realizations of a RRG. Figure [5] analyzes the 
fluctuations about the steady state densities observed 
in an individual run. More precisely, we compute the 
PSNFs in 1.5 x 10 3 parallel long simulation runs on 
RRGs numerically, and then we average. In both figures, 
the same parameter values in the endemic phase are 
used in the corresponding panels to ease a comparison 
of the results. Note that on a finite graph, the only true 
steady state of the SIRS process corresponds to all nodes 
being in susceptible state, that is why the steady states 
of the PA model are compared with the quasi-stationary 
states of large but finite systems. 
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FIG. 4: (Color online) Densities of the infectives and of the 
susceptibles in the PA [gray (green) lines] and averaged nu- 
merical time series (black lines) obtained from Monte Carlo 
simulations of the SIRS stochastic model on a RRG-4 with 
N = 10 6 nodes. All plots were obtained for S = 1, A = 2.5 
and k = 4. Parameters: (a) and (b) 7 = 0.09; (c) and (d) 
7 = 0.04. 

Plots in Fig. 2] show susceptible and infective densities 
as function of time. Numerical solutions of the PA de- 
terministic equations given in the Appendix are plotted 
in gray (green). Black lines are sample averages of the 
densities as obtained from 10 3 realizations of a RRG-4 
with N = 10 6 nodes. For the SIRS model the agreement 
of the solutions of the PA deterministic equations with 
the averaged dynamics on RRGs is sensitive to the rate 
of immunity waning 7, viz it deteriorates with decreasing 



7 . The upper panels in Fig. [4] are plotted for a 
set of parameter values for which the solutions of the 
PA equations reproduce the behavior of the averaged 
times series on a RRG-4 with a good accuracy both in 
the transient and in the quasi-stationary regime. The 
lower panels in Fig. Q] plotted for a smaller value of 7 
with the other parameters being fixed reflect the striking 
difference in the dynamics of the simulations (black 
lines) and of the PA deterministic model [gray (green) 
lines]. Although the steady state values of the densities 
are close, during the transient period the averaged time 
series attenuate more rapidly in the numerical simu- 
lations than the damped oscillations predicted by the PA. 
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FIG. 5: (Color online) Analytical PSNFs of the infectives and 
of the susceptibles in the PA [gray (green) lines] and averaged 
numerical PSNFs (black lines) calculated from Monte Carlo 
simulations of the SIRS epidemic process on a RRG-4 with 
N = 10 6 nodes. All plots were obtained for 5 = 1, A = 2.5 
and k = 4. Parameters: (a) and (b) 7 = 0.09, lin-lin plot; (c) 
and (d) 7 = 0.04, lin-log plot. 

Plots in Fig. [5] are the approximate analytical PSNFs 
[gray (green) lines] given by formula (fT6|) and the av- 
eraged numerical PSNFs (black lines) calculated from 
1.5 x 10 3 replicate simulation runs of the SIRS process on 
a RRG-4 with N = 10 6 nodes. The upper panels in Fig. 
[5] plotted for parameter values used in the upper panels 
in Fig. [4] show the agreement between the analytical and 
numerical PSNFs of the susceptibles and infectives. The 
analytical PSNFs [Eq. (fl"6|)] deduced from the detailed 
stochastic PA model with the transition rates calculated 
on the basis of Eqs. and §5§ approximate the averaged 
numerical PSNFs well in those regions where the PA de- 
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terministic model predicts the same behavior as that of 
the times series obtained from the simulations averaged 
over many realizations. In this region, the typical PSNF 
is a bell-shaped curve indicating that in an individual 
simulation run susceptible and infective densities oscil- 
late in time with frequencies close to the principal fre- 
quency demarked by the maximum of the curve. In fact, 
the large oscillations of the densities are due to a kind 
of internal resonance previously studied for a non-spatial 
predator-prey model [2^ ]. In the lower panels in Fig. 
[5] the amplitude of both the analytical and the numeri- 
cal PSNFs gets much higher, indicating that for smaller 
values of 7 the stochastic fluctuations become more pro- 
nounced and better structured. The fact that the nu- 
merical PSNFs lie significally below those predicted ana- 
lytically is closely related to the stability of the endemic 
equilibrium of the PA equations. As 7 decreases the equi- 
librium's stability gets weaker until it is lost on a critical 
line corresponding to a supercritical Andronov-Hopf bi- 
furcation [271 ] . Accordingly, the approximate analytical 
PSNFs become more and more enhanced until they fi- 
nally diverge on the critical line. Note that in the lower 
panels in Fig. [4] plotted for the same parameters the be- 
haviors of the analytical and numerical global densities 
do not agree too. 

A good correspondence and subsequent divergence be- 
tween the analytical and numerical PSNFs can be under- 
stood from the van Kampen's large- N expansion about 
the endemic equilibrium performed above. The analyti- 
cal expression for the PSNFs, see formula (fT6l) . is a func- 
tion of a; and constant matrices A and B whose elementes 
are expressed in terms of the basic parameters of the 
model A, 7, 5 and k. In the general theory developed 
in Ref. [31|, the matrices can depend on time through 
the node [P(S), P(I)] and pair [P{SI), P(SR), P(RI)} 
densities. However, as we are interested in the analysis 
of the fluctuations in the endemic equilibrium we have 
to substitute for the node and pair variables their sta- 
tionary values [P(S), P(I), P(ST), P{SR), P{RI)} de- 
pending on the parameters. This substitution results in 
time-independent coefficient matrices A and B. It fol- 
lows that the agreement between the analytical and nu- 
merical PSNFs is expectable in the regions where the 
PA deterministic model approximates the transient and 
quasi-stationary global dynamics on RRGs quite well. 
We emphasize that although the coarse grained stochas- 
tic model considered in Ref. [53] exhibits the same qual- 
itative behavior as the detailed model described here, 
the quantitative agreement between the corresponding 
PSNFs is not achieved in that case. Therefore, the de- 
tailed microscopic description is needed to approximate 
the exact stochastic dynamics on RRGs. Moreover, tak- 
ing into account the above analysis it now becomes clear 
that resorting to higher order cluster approximations is 
necessary to describe the behavior of the SIRS stochastic 
model on RRGs for small 7. 



III. DETERMINISTIC AND STOCHASTIC 
FRAMEWORKS BEYOND THE PAIR 
APPROXIMATION 



As mentioned in the previous section, in Ref. 3J] we 
have compared the data of the SIRS stochastic process 
obtained from Monte Carlo simulations on RRGs with 
the solutions of the standard PA equations and have 
shown that the PA describes correctly the global behavior 
of the model in the limit where rate of immunity waning 
7 ^> 1 but it fails to capture the dynamics for 7 <C 1. 
The question is then whether a cluster approximation of 
the next order can explain the suppression of global oscil- 
lations predicted by the PA for 7 <§C 1 and, in particular, 
whether it can describe stationary states correctly. A 
reasonable description of the endemic equilibria as well 
as of the phase diagram of the stationary state calculated 
from numerical simulations is obtained by extending the 
generalized cluster approximation procedure to combina- 
tions of three neighboring nodes or triplets for short. 

In this section, we address this question and more gen- 
erally the problem of the construction of cluster approxi- 
mations of the order q higher than two, that is higher 
than the PA. As a matter of fact, for the SIRS pro- 
cess the time evolution of the g-node joint probabilities 
is governed by a set of first order differential equations 
expressing their time derivatives as linear combinations 
of q- and (q + l)-node joint probabilities. This is due 
to the infection process involving two nodes, susceptible 
and infected, simultaneously. In order to proceed the 
set of equations must be closed. In the standard cluster 
approximation the (g+ l)-node joint probabilities are ra- 
tional functions of the joint probabilities of smaller clus- 
ters of neighboring nodes, appropriately normalized, and 
thus the full set of first order differential equations can 
be obtained. 

Within this perspective the PA is a standard clus- 
ter approximation for q = 2. The TA is obtained for 
q = 3 in a straightforward way by keeping node, pair 
and triplet probabilities as independent variables and ex- 
pressing quadruplet probabilities in terms of them. We 
remind that each equation for the probability evolution 
of a particular cluster configuration is derived by con- 
sidering all transitions leaving or entering it. Using the 
notation of the Appendix, in the TA of the SIRS process 
the equation for, for example, P(RI)(t) reads as 



dP(RI) 
dt 



5P(II) -(S + j)P(RI) + X(k - l)P(RSI) 



The first term on the right-hand side is due to the tran- 
sition of 77 pairs to RI pairs occurring at rate S. The 
second [third] term corresponds to recovery [loss of im- 
munity] of I [R] node in RI pairs that transit to RR 
[SI] pairs at rate 6 [7] . Finally, a triplet RSI changes to 
RII with rate A such that a pair RS is changed to RI. 
Since S node of RS pair has (k — 1) free neighbors that 
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can infect it, this factor is included in the fourth term 
of the equation. Keeping again the same node [P(S), 
P(I)} and pair probabilities [P{SI), P(SR), P(RI)} as 
independent variables the equation can be written in the 
following form: 



dP(RI) 
Jt 



S[P(I) P(SI)} - 
A(fc - l)P(RSI) 



{25 + 1 )P{RI) 



The equations for triplet probabilities are obtained in 
a similar way. In general, as far as clusters of more 
than two neighboring nodes are concerned, these split 
into two distinct classes, open and closed, and both have 
to be taken into account by considering the probabilities 
of finding an open and a closed configuration separately. 
By definition, an open cluster of a given size does not 
contain any loop while a closed cluster, on the contrary, 
necessarily contains at least one loop. Thus, for exam- 
ple, a triplet cluster can be closed forming a triangle or 
open forming a linear chain of three nodes. However, 
the number of short loops is small for large RRGs with 
small k [5^1 . For instance, the analytical estimates for the 
mean Ni and variance Var(Ni) of the number of loops of 
length I = 3, 4 in a RRG-4 in the thermodynamic limit 
are N 3 = Var(N 3 ) = 4.5, N 4 = Var(N 4 ) = 10.125 HH- 
The RRGs with k — 3 are even more sparse, so that the 
number of short loops is even smaller. Based on these 
results we neglect the presence of small closed clusters in 
RRGs. Thus in the TA of the SIRS process, we comple- 
ment Eq. (|27|) . see the Appendix, in which the probabil- 
ities of the triplets will be retained, with the equations 
for the triplet probabilities considering only open clusters 
of neighboring nodes. For instance, using the transition 
rules of the SIRS process the evolution equation for the 
probability P(RRI) of finding a random triplet in state 
RRI becomes: 



dP{RRI) 
Jt 



5[P{RII) + P{IRI) - P(RRI)} 
2 1 P{RRI) + \{k - l)P(RRSI) 



where P(RRSI) is the probability of an open linear 
quadruplet in state RRSI. The deduction of other triplet 
equations is straightforward. 

To write the final closed set of differential equations in 
the TA we have: (a) to choose independent triplet vari- 
ables whose number is reduced due to basic conservation 
relations; (b) to approximate quadruplet probabilities in 
terms of the node, pair and triplet ones. 

We discuss these two questions by order. With the use 
of probabilistic relations 



P(a(3) = Y, P (^X) = Y P (* a ® 



and reflection symmetries 



P(ap X ) = P(xM 



(17) 



(18) 



b) 



a — 0-x-Z 



a — P — X 
I 



FIG. 6: (a) Linear and (b) "T-like" quadruplets both belong 
to the class of open quadruplet configurations, i.e. clusters of 
four nodes that do not contain any loop. 



where a,(3,x € {5, /, i?}, the number of independent 
triplet variables yields 9 and thus the total number of 
the TA equations is 14. Note that in finite RRGs the 
probability of triplets at given time equals: 



P(ap X ) 



k{k- l)N 
k(k- 1)N 



if a ^ X 
if a = x 



(19) 



where we used the obvious notation for triplets. In this 
case, the factor of 2 is due to double counting of af3a 
triplets. From Eqs. (fl~7|) -(|19 p . the constraints on the 
number of pairs and triplets of the following type can be 
found [compare with Eqs. (fl]) and @]: 



(k - l)m a p 



m a /3 x + 2m a p c 



(20) 



A closer inspection of the equations for the triplets 
shows that probabilities of two types of open quadruplet 
configurations occur: linear quadruplets as in Fig. [6] (a) 
and "T-like" quadruplets as in Fig. [6] (b). Note that 
open triplet configurations can only be linear. With the 
assumption of uncorrelated triplets, linear quadruplets 
are approximated as follows. The joint probability of a 
quadruplet in state a(3x£, equals the joint probability to 
find a triplet in state afix multiplied by the conditional 
probability to have a node in state £ neighboring with 
X/3a triplet: 



P(aPxO = P(apx)Q(Z\xPa) • 



(21) 



Neglecting the influence of a node on £ node, the distri- 
bution of £ nodes in the neighborhood of xP a triplets is 
approximately equal to that in the neighborhood of xP 
pairs: 



Qitlxfa) « Q(£\x0) 



(22) 



Substitution of Eq. (|2"2"|) into Eq. (|2"Tj) yields the closure 
assumption for linear quadruplets in the standard TA: 



P(apxO 



P(aPx)P(PxO 

P(Px) 



(23) 
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Note the resemblance of Eq. ([23)1 with the closure as- 
sumption for open triplets in the standard PA: 



P(aPx) 



P(a/3)P(J3 X ) 



(24) 



In both formulas the product of the probabilities (of pairs 
in the PA and triplets in the TA) is divided by the prob- 
ability of the configuration in which they overlap (a node 
and a pair, respectively). Reasoning in the same way, 
a closure assumption for a "T-like" cluster depicted in 
Fig. [5] (b), can be obtained. For example, considering 
a pair in state /3£ as the overlapping configuration the 
probability a "T-like" quadruplet becomes: 



P 



£ ) 



P(a0Z)P(Z0 X ) 



(25) 



Clearly, this closure assumption is not unique as any of 
the three pairs /3a, /3%, and /3£ joining in the node (3 is 
suitable as an overlapping configuration. However, they 
all coincide in the limit when triplets are considered to 
be formed by uncorrelated pairs: 



( a f X ) 



P(0a)P(l3x)P(PO 
PW 



(26) 



Notwithstanding the closure assumption for "T-like" 
quadruplets, Eq. ([2"5| . is not unique it includes more 
information since the probability of a quadrupet is ex- 
pessed in terms of the probabilities of triplets and pairs, 
unlike Eq. (|26p where the same probability depends on 
the probabilities of pairs and nodes. There is no a priori 
reason that the TA model with Eq. ([23)) should give a 
better approximation to Monte Carlo simulations than 
the same model with Eq. (|26[) and direct numerical cal- 
culation of the distribution of nodes, pairs, triplets and 
quadruplets in the simulations is required to check this 
point. However, indirect evidence of this comes from the 
comparison of the TA model with both types of closure 
assumptions for "T-like" quadruplets and the stochastic 
SIRS process on RRGs in the regions where the PA model 
fails. We have found that the closure assumption given 
by Eq. ([21))) gives a quantitative improvement over the 
closure given by Eq. ([2"6")) both for the stationary and for 
the time-dependent behaviors when it is faced with the 
exact stochastic dynamics on RRGs. Qualitatively, the 
behavior of both TA models is similar. For k = 4 no sta- 
ble oscillatory behavior is observed for small 7, and for 
k = 3 the oscillatory phase becomes much smaller than in 
the PA model. The only difference we have been able to 
identify without performing the full linear analysis and 
calculating the whole phase diagrams of the TA deter- 
ministic models is that in the TA model with Eq. ([25)) 
the oscillations are suppressed faster than in the same 
model with Eq. ([26)) . making the former model a better 
approximation for the global behavior observed in the 
simulations on RRGs. 



In Figs. [7] and [8] we compare the three models of 
the SIRS process, namely the TA-SIRS model with the 
closure assumptions for linear and "T-like" quadruplets 
given by Eqs. ([2"3")) and ([21))) [solid gray (green) lines], 
the PA- SIRS model given in the Appendix [dashed black 
(blue) lines] and the results of stochastic simulations 
(solid black lines) on RRGs with k = 3 and k = 4, re- 
spectively. In the plots where the solid black lines can- 
not be distinguished they are superimposed by the solid 
gray (green) lines. The differential equations were in- 
tegrated numerically using the 4th order Runge-Kutta 
algorithm, and for each set of initial conditions and pa- 
rameters the simulations were averaged over an ensemble 
of 10 3 RRGs of given degree with N = 10 6 nodes. The 
properties of the fluctuations around this averaged be- 
havior are described by the numerical PSNFs. For all pa- 
rameter values considered below these numerical PSNFs 
are resonant-like as in Fig. \5\ 
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FIG. 7: (Color online) Time evolution of susceptible (left pan- 
els) and infective (right panels) densities for parameter values 
in the endemic region of the phase diagram of the PA model as 
predicted by the TA model [solid gray (green) lines], the PA 
model [dashed black (blue) lines] and the results of stochas- 
tic simulations (solid black lines) on a RRG-3 with N = 10 6 
nodes. All plots were obtained for 5=1, A = 15 and k = 3. 
Parameters: (a) and (b) 7 = 2; (c) and (d) 7 = 0.2; (e) and 
(f) 7 = 0.05. 



In Fig. [7] the evolution of susceptible (left panels) and 
infective (right panels) densities as a function of time 
is shown for three sets of parameter values that corre- 
spond to constant infection rate A = 15 and decreasing 
rate of immunity waning 7 = 2, 0.2, 0.05 (from top to 
bottom) . The values of the parameters are chosen so as 
to reflect the behavior of the models in three different 
phases in the endemic region of the phase diagram of the 
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PA model 34( . As regards the asymptotic behavior of the 
PA-SIRS equations, these phases are associated with the 
asymptotically stable nodes, asymptotically stable foci 
and limit cycles (from top to bottom). In the upper 
panels for the typical value 7 = 2 we used, the steady 
state of the system is a stable node both for the PA and 
for the TA model. The solutions of both deterministic 
models are almost coincident with the averaged densities 
obtained from the stochastic simulations for this whole 
region. The agreement between the deterministic models 
deteriorates in the region where both predict a behavior 
corresponding to a stable focus as can be seen in the mid- 
dle panels where we used 7 = 0.2. The steady state of 
the PA equations is different from the quasi-stationary 
state observed in the simulations while the TA model 
describes the dynamics accurately both in the transient 
and in the steady state regime. The bottom panels with 
7 = 0.05 illustrate that in the region where the PA model 
exhibits stable oscillatory behavior both the TA solutions 
and the simulation trajectories show damped oscillations 
towards a nontrivial equilibrium. The steady state den- 
sities given by the TA equations and the quasi-stationary 
densities calculated from the stochastic simulations on a 
RRG-3 are equal. However, in the transient regime we 
observe that the trajectories approach the steady state in 
a slightly different manner demonstrating a higher damp- 
ing in the stochastic simulations. Also, for k — 3 the os- 
cillatory phase still persists in a very small region in the 
endemic phase in the TA but once more this result is not 
confirmed by the simulations for the same parameter val- 
ues [results not shown]. Sustained oscillations, instead of 
resonant fluctuations, would show up as multiple peaks 
in the numerical PSNFs, which are not observed. This 
suggests that a complete description of the global behav- 
ior of the SIRS on a RRG-3 requires even more elaborate 
approximations and that the global oscillations predicted 
by both the PA-SIRS and the TA-SIRS in the thermo- 
dynamic limit are an artifact of the models. 

Fig. [5] illustrates the data as in Fig. [7] with A = 2.5, 
and 7 = 2.5,0.1,0.025 (from top to bottom) for k = 4 
and RRG-4. These values are chosen as before to rep- 
resent the different phases of the PA diagram associated 
with the asymptotically stable nodes, asymptotically sta- 
ble foci and limit cycles (from top to bottom). We ob- 
serve the same comparative behavior of the models as 7 
decreases. The only difference is that in the case of k = 4 
the oscillatory behavior in the TA has not been identified. 
Everywhere in the region of stable oscillations predicted 
by the PA model, the TA model reaches the steady state 
through damped oscillations that are associated with the 
existence of a stable fixed point, namely a stable focus. 
Such a dynamics of the TA-SIRS model is confirmed by 
the results of the simulations, see the bottom panels in 
Fig. E 

Regarding the fluctuations, the same procedure that 
was used in the previous section to derive an approximate 
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FIG. 8: (Color online) The same data as in Fig. [TJwith 5 = 1, 
A = 2.5 for k = 4 and RRG-4 with N = 10 6 nodes. Param- 
eters: (a) and (b) 7 = 2.5; (c) and (d) 7 = 0.1; (e) and (f) 
7 = 0.025. 



analytical expression for the PSNFs based on the PA can 
be extended to the TA in the region where the PA model 
fails. An explicit construction of a detailed stochastic TA 
model is cumbersome but straightforward, and once the 
master equation is written down the analytical PSNFs 
can be computed using van Kampen's system size ex- 
pansion and then the general formula (|16p. According 
to the results for the PA model, we expect the analyti- 
cal PSNFs obtained in this way to approximate well the 
numerical power spectra in the parameter regions of the 
middle panels in Fig. [7J We also expect a fair approxi- 
mation in the parameter regions of the bottom panels of 
the same figure. 



IV. DISCUSSION AND CONCLUSIONS 

The standard PA is known to perform poorly for 
lattice-based stochastic models but it gives in general 
good results for large RRGs due to the local tree-like 
structure of these graphs. Taking a simple epidemic 
model as an example, we have used the PA to derive 
the master equation of the corresponding stochastic pro- 
cess on a RRG and an approximate analytical expression 
for the power spectrum of the fluctuations in the quasi- 
stationary state. 

We have checked the agreement of the analytical power 
spectrum in the PA against numerical simulations and 
found that whenever the behavior of the system in the 
thermodynamic limit is well described by the PA deter- 
ministic equations, the analytical power spectrum also 
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describes accurately the fluctuations observed in long 
simulations. 

This happens in a large region of parameter space. 
However, as 7 approaches the phase boundary 7 = 
the quality of the PA deteriorates, and it is necessary 
to switch to higher order cluster approximations in or- 
der to obtain even qualitative agreement between the 
model equations and the simulations. We have shown 
that a TA with a standard closure assumption yields an 
accurate description of the behavior of the system in the 
thermodynamic limit in a 7 range where the PA breaks 
down. For finite systems, the fluctuation power spectrum 
can be computed analytically as before from the master 
equation of the stochastic process that corresponds to the 
TA. 

For small values of 7, long simulations require very 
large system sizes. For the smallest 7 we have explored 
we found indications of the breakdown of the TA, and 
that clusters of order higher than three would have to 
be considered. As the order of the cluster approximation 
increases, however, the 'no loop' assumption that is an 
ingredient of the construction of the detailed stochastic 
model becomes less accurate. Extension of this method 
to smaller values of 7 through higher order clusters would 
have to be combined with more complicated closure as- 
sumptions. 

The parameter range explored in this paper is rele- 
vant to childhood infectious diseases modeling. Pub- 
lished estimates for the epidemiological parameters of 
measles, whooping cough, rubella and chicken pox, see 
Ref. Q , correspond in the SIRS model to 7 in the range 
0.003 ~ 0.02 assuming that the average immunity wan- 
ing period can be taken as the typical duration of basic 
school. Our results show that the oscillatory phase that 
implicitly spatial models such as the PA-SIRS and the 
TA-SIRS exhibit in the thermodynamic limit cannot be 
directly related with the recurrent epidemic peaks found 
in many data sets [H, [2|| . However, they also show that 
once stochastic effects are taken into account, the model 
predicts a well defined bell-shaped high amplitude fluc- 
tuation spectrum reproducing the qualitative features of 
typical time patterns of real data for this class of endemic 
diseases. 
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APPENDIX 



Throughout the main text the node and pair proba- 
bilities are denoted as P(a) and P(a(3), where the small 
Greek letters stand for states S, I, and R. We extend 
this notation for open linear triplet and quadruplet prob- 
abilities, P(a/3x) and P(af3x£), respectively, while in the 
probabilities of open " T-like" quadruplets the clusters are 
depicted explicitly 

Substituting P(a/?x) = P{aP)P((3x)/ 'P(P) in the set 
below yields the PA-SIRS deterministic equations [271 ]: 



dP(S) 

dt 
dP(I) 
dt 
dP{SI) 

dt 



dP(SR) 
dt 

dP{RI) 
~dt 



= 7 [1 - P{S) - P(I)} - kXP(SI) , (27) 

= kXP(SI) - 5P{I) , 

= jP(RI) - SP(SI) - XP(SI) + 

+ \(k - 1)[2P(SSI) + P(RSI) - P(SI)} , 

= SP(SI) - X(k - 1)P(RSI) + 

+ 7 [1 - P(S) - P(I) - P(RI) - 2P(SR)} , 

= 6 [P(I) - P(SI) - 2P(RI)} - jP(RI) + 

+ \(k - 1)P{RSI) . 



The steady state solutions of the PA-SIRS equations 
can be obtained analytically. Let P(S), P(I) and P(SI), 
P(SR), P(RI) denote the endemic steady state values of 
the node and pair densities of the PA-SIRS model, then 
the Jacobian matrix A of the linearized system is written 
as 



-7 


-7 


-kX 





-5 


kX 


C0C2 





C 3 



A = 



-7 + Co 
-Co 








Pisp 



P(SR) 
-27- 1 



(28) 



P(SR) 



P(SI) P(SR) 



- 7 -25 
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where we introduced the constants 

(k - 1)XP(SI)P(SR) 



Co 



C x = 



(k - 1)XP(SI)P(SR) 



P(S) 



C-2 



, Hsi) 

'P{SR) 



and 



c ^^ k - 2)x - 5 -^{pW) + Wm) • 

In the TA of the SIRS process, we complement Eq. ([77]) 
in which the probabilities of the triplets are retained with 
the equations for the triplet probabilities as described in 
the main text. 
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